Using the sample information from the phase I trial, we can calculate power for a phase II trial.
We wish to identify sample sizes for high dose and standard dose (in a 2:1 ratio) that will be able to detect a positive difference in the probability $p$ of a $\ge 4$ fold titer rise to H3N2 influenza. That is, we wish to be able to detect:
$$\delta = p^{(HD)} - p^{(SD)} > 0$$with 80% power. This analysis will be conducted using the proper Bayesian approach (Spiegelhalter et al. 2004). The advantage of this approach for power analysis over a frequentist power analysis is that it does not assume that the underlying probabilities are known. Failure to do so will tend to underestimate variance and result in underpowered experiments.
A sample leading to a successful trial (one in which a difference is detected) under a particular sample size $n$ is:
$$E_n = \left\{\mathbf{x}; Pr( \delta > 0 | \mathbf{x} ) \ge 0.95\right\}$$so the probability of a successful trial is:
$$\psi(n) = P_{X^{(n)}}(X^{(n)} \in E_n)$$which is the expected Bayesian power.
In our case, the power is:
$$\psi(n) = \Phi\left(\frac{\pi_{hd} - \pi_{sd}}{\sqrt{\pi_{hd}(1-\pi_{hd})/n_{hd} - \pi_{sd}(1-\pi_{sd})/n_{sd}}}\right)$$where the $\pi_i$ are the posterior probabilities, parameterized by a beta distribution with a binomial likelihood for the data.
We will use the sample data from the phase I trial (GiaQuinta et al. 2014) to construct prior distributions for $p^{(HD)}$ and $p^{(SD)}$.
%matplotlib inline
import numpy as np
import seaborn as sns
import pymc3 as pm
import theano.tensor as tt
import matplotlib.pyplot as plt
from theano.tensor.shared_randomstreams import RandomStreams
xvals = np.linspace(0, 1)
hd_dist = pm.Beta.dist(alpha=13, beta=11).logp(xvals)
sd_dist = pm.Beta.dist(alpha=3, beta=14).logp(xvals)
plt.plot(xvals, np.exp(hd_dist.eval()))
plt.plot(xvals, np.exp(sd_dist.eval()))
[<matplotlib.lines.Line2D at 0x78ee8af53240>]
# Standard normal from complimentary error function
Φ = lambda x: 0.5*tt.erfc(-x/tt.sqrt(2))
def power_analysis(n):
# Allocate samples 2:1
n_hd = int(2*n/3)
n_sd = n - n_hd
with pm.Model() as model:
# probabilities of >= 4-fold-rise taken from Beta priors constructed from phase I data
π_hd = pm.Beta('π_hd', 12, 10)
π_sd = pm.Beta('π_sd', 2, 13)
# Draw binomial samples from probabilities
rng = RandomStreams(seed=20090425)
#x_hd = pm.Binomial('x_hd', n_hd, π_hd)
x_hd = rng.binomial(n=n_hd, p=π_hd)
#x_sd = pm.Binomial('x_sd', n_sd, π_sd)
x_sd = rng.binomial(n=n_sd, p=π_sd)
# Empirical probabilities
π_hd_post = x_hd/n_hd
π_sd_post = x_sd/n_sd
# Convert difference in probabilities to standard normal
d = pm.Deterministic('d', (π_hd_post - π_sd_post)
/ tt.sqrt(π_hd_post*(1-π_hd_post)/n_hd + π_sd_post*(1-π_sd_post)/n_sd))
# P-value
p = 1 - Φ(d)
# Successful detection
power = pm.Deterministic('power', (p<0.05).astype('int32'))
return model
Use model to estimate power for a range of sample sizes.
n=60:
iterations = 1000
tune = 2000
mod60 = power_analysis(60)
with mod60:
tr = pm.sample(iterations, tune=tune, cores=2)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [π_sd, π_hd] Sampling 2 chains: 100%|██████████| 6000/6000 [00:04<00:00, 1284.20draws/s]
pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]
mean | mc_error | |
---|---|---|
power | 0.891 | 0.006796 |
n=48:
mod48 = power_analysis(48)
with mod48:
tr = pm.sample(iterations, tune=tune, cores=2)
pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [π_sd, π_hd] Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1148.71draws/s] The acceptance probability does not match the target. It is 0.8859333589363031, but should be close to 0.8. Try to increase the number of tuning steps.
mean | mc_error | |
---|---|---|
power | 0.8605 | 0.007625 |
n=42:
mod42 = power_analysis(42)
with mod42:
tr = pm.sample(iterations, tune=tune, cores=2)
pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [π_sd, π_hd] Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1085.89draws/s]
mean | mc_error | |
---|---|---|
power | 0.8485 | 0.008528 |
n=36:
mod36 = power_analysis(36)
with mod36:
tr = pm.sample(iterations, tune=tune, cores=2)
pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [π_sd, π_hd] Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1162.65draws/s]
mean | mc_error | |
---|---|---|
power | 0.834 | 0.007513 |
n=27:
mod27 = power_analysis(27)
with mod27:
tr = pm.sample(iterations, tune=tune, cores=2)
pm.summary(tr, varnames=['power'])[['mean', 'mc_error']]
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [π_sd, π_hd] Sampling 2 chains: 100%|██████████| 6000/6000 [00:04<00:00, 1244.64draws/s]
mean | mc_error | |
---|---|---|
power | 0.7505 | 0.009151 |
The power analysis suggests that a sample size in the low- to mid-30s would be sufficient to detect a difference between high and standard dose with 80% power in a similar population. A sample size over 60 would be required for 90% power.